Thermal quenches in spin ice 
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We study the diffusion annihilation process which occurs when spin ice is quenched from a high temperature 
paramagnetic phase deep into the spin ice regime, where the excitations - magnetic monopoles - are sparse. 
We find that due to the Coulomb interaction between the monopoles, a dynamical arrest occurs, in which non- 
universal lattice-scale constraints impede the complete decay of charge fluctuations. This phenomenon is outside 
the reach of conventional mean-field theory for a two-component Coulomb liquid. We identify the relevant 
timescales for the dynamical arrest and propose an experiment for detecting monopoles and their dynamics in 
spin ice based on this non-equilibrium phenomenon. 



Introduction - There is intense current interest in the study 
of strongly correlated systems hosting fractionalised excita- 
tions, in fields as diverse as magnetism, quantum Hall physics 
and quantum computing, or even the study of (topological) 
band insulators. Such excitations arise against the background 
of highly unusual ground states. 

Recently, we have argued that spin ice - an Ising magnet on 
the pyrochlore lattice - hosts deconfined magnetic monopole 
excitations, which result from the fractionalisation of the high- 
energy local dipole moments 12]. At the moment, the focus of 
theoretical and experimental studies consists of predicting and 
detecting signatures of these excitations J2-01- 

In spin ice, the ground state ensemble is unusual in that it 
exhibits algebraic correlations without representing a conven- 
tional critical point: this Coulomb phase - in the sense of the 
deconfined phase of a U ( 1 ) gauge theory - is a consequence of 
the local constraint that two spins point into each tetrahedron 
and two point out. Indeed, spin ice owes its name to this mag- 
netic version of the Bernal-Fowler ice rules. The Coulomb 
phase is characterised by an emergent gauge field, rather than 
an emergent order parameter; as such, it is a classical example 
of topological order. 

Violating the ice rules by flipping a spin out of a ground- 
state configuration, at a cost in energy of A s f , leads to a pair of 
pointlike defects in the tetrahedra the spin belongs to. These 
two defects are deconfined: they can be separated to an arbi- 
trarily large distance at a finite cost in energy. In the presence 
of long-range dipolar interactions - the model referred to as 
dipolar spin ice below - such defects experience a magnetic 
Coulomb interaction, V(r) = poQ^/^wr), whence the ap- 
pellation magnetic monopoles. Here, po is the vacuum perme- 
ability, and the magnetic charge Q m — 2\j2\/a,d is related to 
the dipole moment of the magnetic ions, \p\, and the distance 
between the centres of adjacent tetrahedra, ad- 

In addition, there is a Coulomb interaction of entropic ori- 
gin, with coupling strength Q\ oc T. It is present even in the 
nearest-neighbour model for spin ice, where the long-range 
dipolar interactions are omitted and Q m = 0. 

Spin ice compounds, such as Ho2Ti2C»7 and Dy2Ti2C<7, are 
thus the first instances of three-dimensional magnets which 
host deconfined fractionalised excitations. The most satisfy- 



ing detection experiment would consist of a direct visualisa- 
tion of a magnetic monopole in bulk spin ice. However, due 
to its small magnetic charge and the fact that single quasipar- 
ticles are hard to come by in bulk systems - even in quantum 
Hall physics, a single fractionalised charge has never been im- 
aged - this has so far proven beyond reach. In Ref. Q1 we 
have shown that a thermodynamic signature of the magnetic 
Coulomb interaction of the monopoles is the presence of a 
liquid-gas transition in a magnetic-field applied in the [111] 
direction, which had already been experimentally observed. 

In this publication, we study the evolution of the monopole 
density after a thermal quench. The description of such non- 
equilibrium dynamics is a worthwhile enterprise in itself, as 
thus far there have been no instances of three-dimensional 
magnets with pointlike elementary excitations, and hence lit- 
tle motivation for their study. However, there has been work 
on the quench dynamics of Coulomb liquids 10, 0], which pro- 
vides the starting point for our analysis. 

Our central result consists of the demonstration that the 
time-dependence of the monopole density after a quench pro- 
vides a distinct signature of not only their pointlike nature but 
also of their magnetic Coulomb interaction. For the nearest- 
neighbour model, we show that mean-field theory applies. We 
analytically account for the simulated time dependence of the 
density without free parameters. In dipolar spin ice, monopole 
bound states appear which can only be annihilated over an en- 
ergy barrier. This leads to a dynamical arrest at low temper- 
atures. This is again borne out by Monte Carlo simulations, 
where the fundamental dynamical move consists of a single 
spin flip as appropriate for the large Ising spins in spin ice |2]. 

The outline of the paper is as follows. We set the stage 
by briefly summarising the annihilation-diffusion physics in 
Coulomb liquids. We address the new features present in 
nearest-neighbor spin ice, before presenting our results on the 
dipolar system. We close with remarks on equilibration in 
spin ice, and how the freezing of bound pairs could be used as 
an experimental technique to achieve measurable monopole 
densities at very low temperatures. 

Diffusion-annihilation in Coulomb liquids - Consider a 
density n±(r) of positive and negative monopoles. As op- 
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positely charged pairs can annihilate, their density obeys 



dn + (r) dn-(r) 



dt 



dt 



-/Cn + (r)n_ (r), 



(1) 



where K. is an appropriate rate constant. In addition, the 
monopoles move deterministically in response to their mu- 
tual forces, and they are subject to diffusion in the presence of 
density inhomogeneities. 

Neglecting density fluctuations, one obtains the mean-field 
solution 



P(t) 



n + + n. 



Po 



1 + ICpot 



for a quench to T = 0, where po 



is the initial density: 



the characteristic timescale for the decay is t/c ~ <2q//C. A 
dynamical bottleneck can arise if there are spatial fluctuations 
in the relative density of positive and negative monopoles 
a(r) — [n+(r) — ra_(r)]/2 (which is unaffected by the sym- 
metric annihilation process) that are not smoothed fast enough 
by the motion of the monopoles. The relevant timescale for a 
particle to move a distance a is given by tq ~ ao/(p>E) ~ 
do/(pq), where p is the monopole mobility and E ~ q/a 2 , is 
the typical strength of the Coulomb field. 

Nearest-neighbor spin ice - This system presents a number 
of special features with respect to ordinary Coulomb liquids, 
which are intricately linked to the existence of the monopoles 
against a backdrop of spin configurations in spin ice. The first 
is a constraint on the possible values of a which follows from 
the fact that the charge density encodes the change in the mag- 
netisation of the sample. The boundedness of the magnetisa- 
tion implies that a cube of volume L 3 can at most accommo- 
date a net charge a ~ L 2 . Thus the long- wavelength Fourier 
components are suppressed as a(q) ~ q 2 . 

The other crucial feature is that the interaction between the 
defects is of a purely entropic nature, due to the weighting of 
the monopole states by the number of spin configurations they 
are compatible with. This interaction has a Coulombic form 



V s (r) = k B T 



r/a,d 



(3) 



1511 . Whereas the 
-)• 0, the mobil- 



where Q 2 ~ 0.35 ± 0.01 can e.g. be obtained from the proba 
bility distribution of the separation of a lone pair of monopoles 
in equilibrium Monte Carlo simulations J^, 
strength of this interaction vanishes as T 
ity p ~ (Qma 2 ,) / (QrhsT) arising from the single spin-flip 
Metropolis dynamics diverges, resulting in a regular T — > 
limit of tq. Here r is the basic unit of time, e.g., the inverse 
of the flip rate of an isolated spin in Monte Carlo simulations. 
In spin ice materials, AC susceptibility measurements seem 
to indicate that r ~ 1 ms II 1 0TI . A simple estimate yields 
JC/a 3 ~ 2q/t, where j e [|, A] : the probability of finding 
two oppositely charged monopoles on neighbouring tetrahe- 
dra is proportional to p 2 , and they annihilate in the next step 
(after time r) if flipping the intermediate spin restores the ice 
rules in both tetrahedra; this probability, which depends on the 



spin correlations, is estimated by g (1 — g being the probabil- 
ity that the two monopoles do not annihilate upon flipping the 
intermediate spin, as illustrated in the left panel of Fig. |2j. 

Our numerical simulations of thermal quenches in nearest- 
neighbour spin ice down to zero temperature, are displayed in 
Fig- El The mean field solution shows quantitative agreement 
with the numerics without any fitting parameters: the fact that 
a(q) ~ q 2 , together with the entropic Coulomb interaction, 
effectively suppress fluctuations in the charge density. 
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FIG. 1. (Color Online) - Monopole density evolution in nearest- 
neighbour spin ice (J = 1 K), after a temperature quench from 
T — 10 K down to T = K. We simulates systems of size 
L = 32, 64, 128, and finite size effects are absent at these time 
scales (estimated error bars are smaller than the symbol size). The 
analytical mean-field result Eq. l(2j is shown for g = 3/4 (dashed 
black line) and g = 9/10 (solid black line). 

Dipolar spin ice - The presence of a magnetic Coulomb in- 
teraction in spin ice leads to further features outside the con- 
ventional picture of Coulomb liquids. First of all, a diverg- 
ing mobility is no longer compensated by a vanishing poten- 
tial energy, and tq — > in the zero temperature limit. This 
indicates that the motion of the monopoles does not follow 
linear response but rather the monopoles move along the lo- 
cal field direction at the maximum speed permitted by micro- 
scopic constraints, namely one step in time r. We thus expect 
monopoles to find each other very efficiently, and therefore a 
decay of p which is at least as fast as in the nearest-neighbour 
case. 

This is, however, not what happens. 

The interplay between long-range interactions and con- 
straints imposed by the underlying spin degrees of freedom 
leads to the formation of non-contractible monopole pairs, 
and the system exhibits a dynamical arrest. Indeed, not all 
nearest-neighbour monopole-antimonopole pairs can be an- 
nihilated by flipping the shared spin (see Fig. [2] left panel). 
Annihilation can then take place only if the two monopoles 
separate and meet again elsewhere in the lattice. However, 
due to their magnetic Coulomb interaction, there is an energy 
barrier for such process, leading to an activated Arrhenius be- 
haviour in the monopole density relaxation. 
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FIG. 2. (Color Online) - Example of a non-contractible monopole- 
antimonopole pair (left panel). The shortest path that can lead to their 
annihilation is a hexagonal loop, provided the spins along the path are 
oriented appropriately (right panel). One can see explicitly that the 
two monopoles must separate before they are allowed to annihilate, 
resulting in a Coulomb energy barrier for the process. 



The smallest possible energy barrier determines the long 
time behaviour in the system. This is given by an elemen- 
tary move where the monopoles of a bound pair annihilate 
around one of the adjacent hexagonal loops in the lattice (see 
Fig. 12 right panel). Two of the five spin flips involved in 
such process increase the distance between oppositely charged 
monopoles. A rough estimate for the concomitant energy 
gaps (see Ref. [2) is given by the Coulomb interaction be- 
tween the magnetic charges. From the nearest neighbour 
value E nn = /ioQf n /4:irad — 3.06 K, we obtain the bar- 
rier to hop to second neighbour distance a2„ = y/8/3ad, 
Ai = E nn (l — a,d/a2n) — 1-19 K, and the barrier to hop 
from second to third neighbour distance a^ n = y/ll/Saa, 
A 2 = E nn (a d /a 2n - a2„/a 3 „) = 0.28 K, leading to an over- 
all energy barrier of A = Ai + A 2 = 1.47 K. In practice, the 
energy cost of a spin flip varies due to the effectively random 
fields set up by nearby bound pairs, leading to a broadened 
distribution of the Aj. 

We ran extensive numerical Monte Carlo (MC) simulations 
treating the long range dipolar interaction via the Ewald sum- 
mation technique, [11] and using the Waiting Time Method 
(WTM) jl2ll with single spin flip updates to access the long 
time regime [13]. We prepare the system at equilibrium at the 
initial temperature of 10 K; we then set the temperature to its 
quench value at time t = 0, and we start the measurements. 

The defect density either reaches its equilibrium value very 
quickly, (for T > 0.4 K), or a significant deviation from 
power law decay appears (T < 0.4 K) due to the activated 
behaviour induced by the non-contractible bound pairs, as il- 
lustrated in Fig. [3] 

A (temperature independent) Gaussian distribution of en- 
ergy barriers A peaked around 1.47 K, with a variance 
0.01 K 2 , leads to a probability distribution V(Q) of single 
hexagon decay times 6, and hence a (normalised) defect den- 
sity p(t) = 1- Jo V{&) dQ. The resulting curves p(t) are 
compared with the numerical ones in the inset of Fig. [3] No- 
tice the good agreement over more than 20 orders of mag- 
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FIG. 3. (Color Online) - Numerical simulations of thermal quenches 
in dipolar spin ice (system size L = 8, i.e., 8192 spins, and simula- 
tion parameters for Dy2Ti207 as in Ref. (3). The curves show the 
total density of defects p per tetrahedron as a function of Monte Carlo 
time in units of Monte Carlo steps (1 attempt per spin), for quenches 
from T = 10 K, down to T = 0.025 K (red), T = 0.04 K (blue), 
T = 0.05 K (green), T = 0.075 K (magenta), T = 0.1 K (cyan), 
T = 0.125 K (yellow), T = 0.15 K (black), T = 0.4 K (red), 
T — 0.5 K (blue), and T = 0.6 K (green) - appearing in order from 
right to left. Inset: long time behaviour of p normalised by its plateau 
value ppiatcau, compared to the phenomenological model discussed 
in the text (thin black lines). 



nitude in t, for the different values of the quench tempera- 
ture. Clearly, this phenomenological model captures the fun- 
damental physics underlying the dynamical arrest in thermal 
quenches. 

To further confirm this scenario, we explicitly determined 
the density of monopoles forming non-contractible pairs, as 
well as the density of contractible defect pairs (i.e., pairs 
where flipping the intermediate spin lowers the number of 
defects in the system). The result is illustrated in Fig. [4] 
for a given quench temperature. One can see that the initial 
decay ends when there are essentially no contractible pairs 
left in the system (magenta curve falling below 1/N t , where 
Nt = 8L 3 is the total number of tetrahedra in the lattice). 
From thereon, the total defect density is essentially given by 
monopoles forming non-contractible pairs. 

The defect density decay approaching the plateau is cap- 
tured by a diffusion process where oppositely charged par- 
ticles (A, B) can either annihilate (0) or fuse into a non- 
contractible pair (D) (Fig. [5] left panel). For quenches to 
very low temperatures, the non-contractible pairs can be ap- 
proximated as frozen unless another single particle annihilates 
one member of the pair, thus freeing the other one (Fig. [5] 
right panel). In a simple mean field model, one finds a sur- 
viving population of non-contractible pairs D, provided the 
single particle density decays faster than 1/t, as is the case 
in our simulations. Indeed, the resulting time-dependence 
of the total and non-contractible particle densities is in good 
qualitative agreement with the numerical results illustrated in 
Fig. S] 11511 . On the longest timescales, annihilation of non- 
contractible pairs D — > around hexagonal loops terminates 
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FIG. 4. (Color Online) - Numerical simulation of a thermal quench 
down to T = 0.125 K (system size L = 8). The red curve shows 
the total density of defects per tetrahedron p, while the blue and the 
rapidly decaying magenta curves correspond to the density of defects 
forming non-contractible pairs p nc and contractible pairs p c , respec- 
tively. 
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FIG. 5. (Color Online) - Schematic illustration of the reaction pro- 
cesses in the mean field model used to describe the very low temper- 
ature limit of thermal quenches in spin ice. 

the plateau. 

Equilibration timescales, and experiment - The dominant 
dynamics in spin ice at low temperatures consists of hopping 
monopoles. These correspond to single spin flips which do 
not incur a cost for violating the ice rules. As the temperature 
is lowered to zero, the monopole density vanishes and spin ice 
freezes completely. At finite temperatures the low density of 
monopoles leads to an exponentially large timescale (in fact, 
possibly super-exponentially large B2l ll0ll ) which grows faster 
than the timescale governing the monopole density. Upon 
cooling, there is a "fast" process responsible for the thermali- 
sation of the energy (i.e., the monopole density), and a much 
slower process that equilibrates the spin correlations. We be- 
lieve that this mechanism explains why certain quantities such 
as the energy seem to equilibrate at temperatures where the 
magnetisation has long fallen out of equilibrium. It would be 
on the long time scales of the slower process that any magnetic 
order would be established 11611 . 

In addition, we note that non-contractible bound monopoles 
are practically absent upon heating from an equilibrium state. 
There is thus an asymmetric approach to equilibrium in 
the monopole density at a given temperature, depending on 
whether we use heating or cooling thermal quenches. 

The most elegant way to measure monopole densities 



would be via zero-field NMR on the oxygen nuclei at the cen- 
tre of the tetrahedra, which experience different field strengths 
and corresponding fluctuation rates when the tetrahedra host a 
monopole 111711 . The aim is to have a sufficiently large density 
of monopoles to yield a measurable signal, while preventing 
the monopoles from moving around too quickly, thus spoil- 
ing the measurement. In thermal equilibrium, temperatures 
low enough for the second condition to be satisfied result in 
exceedingly small monopole densities. Our results suggest 
that a temperature quench on time scales sufficiently shorter 
than the time to develop the dynamical arrest plateau in Fig. [3] 
(~ 10-100 ms) could be used to induce a monopole-rich state 
(p > 10 -2 per tetrahedron), where motion at short times is 
obstructed by the formation of non-contractible pairs. 

Finally, the method of choice for imaging spin correlations 
is neutron scattering [3-5]. As the non-contractible monopole 
pairs remain bound on long timescales, the concomitant short- 
range correlations should be visible in the neutron scattering 
cross section. 
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